data_root_water$Bodentiefe <- as.factor(data_root_water$Bodentiefe)
root_water_response <- data_root_water$Prozentualer_Gewichtsverlust_Wurzeln
root_water_parameters <- data_root_water %>%
dplyr::select(where(is.numeric)) %>%
dplyr::select(
-Plot,
-Prozentualer_Gewichtsverlust_Wurzeln
)
root_water_gam_results <- data.frame(
water_variable = names(root_water_parameters),
n_used = NA,
AIC_global = NA,
AIC_depth_smooth = NA,
delta_AIC = NA,
chosen_model = NA,
r_squared_chosen = NA,
p_value_smooth_chosen = NA
)
for (i in seq_along(root_water_parameters)) {
water_values_current <- root_water_parameters[[i]]
root_water_gam_input <- data.frame(
response_root = root_water_response,
water_current = water_values_current,
Bodentiefe = data_root_water$Bodentiefe
)
root_water_gam_input <- root_water_gam_input[complete.cases(root_water_gam_input), ]
root_water_gam_results$n_used[i] <- nrow(root_water_gam_input)
if (nrow(root_water_gam_input) < 12) next
if (length(unique(root_water_gam_input$water_current)) < 5) next
if (length(unique(root_water_gam_input$Bodentiefe)) < 2) next
mod_global <- try(
mgcv::gam(
response_root ~ s(water_current, k = 4) + Bodentiefe,
data = root_water_gam_input,
method = "REML"
),
silent = TRUE
)
if (names(root_water_parameters)[i] == "mean_summer" &&
!inherits(mod_global, "try-error")) {
mod_root_summer_global <- mod_global
}
mod_depth <- try(
mgcv::gam(
response_root ~ s(water_current, k = 4) +
s(water_current, by = Bodentiefe, k = 4) +
Bodentiefe,
data = root_water_gam_input,
method = "REML"
),
silent = TRUE
)
if (inherits(mod_global, "try-error") || inherits(mod_depth, "try-error")) next
if (names(root_water_parameters)[i] == "mean_summer") {
mod_root_summer_global <- mod_global
mod_root_summer_depth  <- mod_depth
plot_data_root_summer  <- root_water_gam_input
}
AIC_global <- AIC(mod_global)
AIC_depth  <- AIC(mod_depth)
root_water_gam_results$AIC_global[i] <- AIC_global
root_water_gam_results$AIC_depth_smooth[i] <- AIC_depth
root_water_gam_results$delta_AIC[i] <- AIC_depth - AIC_global
if (AIC_global <= AIC_depth) {
best_mod <- mod_global
root_water_gam_results$chosen_model[i] <- "global_smoother"
} else {
best_mod <- mod_depth
root_water_gam_results$chosen_model[i] <- "depth_specific_smoother"
}
root_water_gam_results$r_squared_chosen[i] <- summary(best_mod)$r.sq
pvals <- summary(best_mod)$s.table[, "p-value"]
root_water_gam_results$p_value_smooth_chosen[i] <- min(pvals, na.rm = TRUE)
cat("\n========================================\n")
cat("GAM für Water-Level-Variable:", names(root_water_parameters)[i], "\n")
cat("AIC global:", round(AIC_global, 3), "\n")
cat("AIC depth-specific:", round(AIC_depth, 3), "\n")
cat("-> gewählt:", root_water_gam_results$chosen_model[i], "\n")
cat("========================================\n\n")
print(summary(best_mod))
cat("\n")
gam.check(best_mod)
}
root_water_gam_results_sorted <-
root_water_gam_results[order(-root_water_gam_results$r_squared_chosen), ]
root_water_gam_results_sorted
AIC(mod_root_summer_global)
AIC(mod_root_summer_depth)
par(mfrow = c(1, 2))
plot(mod_root_summer_global, shade = TRUE, main = "Globales Modell")
plot(mod_root_summer_depth, shade = TRUE, main = "Tiefenspezifisches Modell")
water_current = seq(
min(plot_data_root_summer$water_current, na.rm = TRUE),
max(plot_data_root_summer$water_current, na.rm = TRUE),
length.out = 200
),
##### BRT ########
root_pairs_table <- basis_data %>%
dplyr::select(
Plot, study_site, Bodentiefe, Pseudoreplikat,
Prozentualer_Gewichtsverlust_Wurzeln
) %>%
dplyr::filter(Pseudoreplikat %in% c("A","B")) %>%
tidyr::pivot_wider(
names_from  = Pseudoreplikat,
values_from = Prozentualer_Gewichtsverlust_Wurzeln,
values_fn   = mean   # falls es Duplikate gibt, werden sie gemittelt
) %>%
dplyr::mutate(
diff_AB     = A - B,
abs_diff_AB = abs(diff_AB),
mean_AB     = rowMeans(cbind(A, B), na.rm = TRUE),
complete_pair = !is.na(A) & !is.na(B)
) %>%
dplyr::arrange(Plot, Bodentiefe)
root_pairs_table
root_pairs_table_complete <- root_pairs_table %>%
dplyr::filter(complete_pair)
root_pairs_table_complete
rhizome_pairs_table <- basis_data %>%
dplyr::select(Plot, study_site, Bodentiefe, Pseudoreplikat, Prozentualer_Gewichtsverlust_Rhizome) %>%
dplyr::filter(Pseudoreplikat %in% c("A","B")) %>%
pivot_wider(names_from = Pseudoreplikat, values_from = Prozentualer_Gewichtsverlust_Rhizome, values_fn = mean) %>%
mutate(diff_AB = A - B, abs_diff_AB = abs(diff_AB), mean_AB = rowMeans(cbind(A, B), na.rm = TRUE),
complete_pair = !is.na(A) & !is.na(B)) %>%
arrange(Plot, Bodentiefe)
prod_pairs_table <- basis_data %>%
dplyr::select(Plot, study_site, Bodentiefe, Pseudoreplikat, BM_accumulation_belowground_per_m2) %>%
dplyr::filter(Pseudoreplikat %in% c("A","B")) %>%
pivot_wider(names_from = Pseudoreplikat, values_from = BM_accumulation_belowground_per_m2, values_fn = mean) %>%
mutate(diff_AB = A - B, abs_diff_AB = abs(diff_AB), mean_AB = rowMeans(cbind(A, B), na.rm = TRUE),
complete_pair = !is.na(A) & !is.na(B)) %>%
arrange(Plot, Bodentiefe)
pfp_pairs_table <- basis_data %>%
dplyr::select(Plot, study_site, Bodentiefe, Pseudoreplikat, peat_formation_potential_depth_wise) %>%
dplyr::filter(Pseudoreplikat %in% c("A","B")) %>%
pivot_wider(names_from = Pseudoreplikat, values_from = peat_formation_potential_depth_wise, values_fn = mean) %>%
mutate(diff_AB = A - B, abs_diff_AB = abs(diff_AB), mean_AB = rowMeans(cbind(A, B), na.rm = TRUE),
complete_pair = !is.na(A) & !is.na(B)) %>%
arrange(Plot, Bodentiefe)
all_pairs_table <- dplyr::bind_rows(
root_pairs_table  %>% dplyr::mutate(response = "roots"),
rhizome_pairs_table %>% dplyr::mutate(response = "rhizomes"),
prod_pairs_table  %>% dplyr::mutate(response = "production"),
pfp_pairs_table   %>% dplyr::mutate(response = "pfp")
) %>%
dplyr::select(response, Plot, study_site, Bodentiefe, A, B, diff_AB, abs_diff_AB, mean_AB, complete_pair) %>%
dplyr::arrange(response, Plot, Bodentiefe)
all_pairs_table
pairs_complete <- all_pairs_table %>%
dplyr::filter(complete_pair)
between_plot_summary <- pairs_complete %>%
dplyr::group_by(response, Bodentiefe) %>%
dplyr::summarise(
n_plots = dplyr::n_distinct(Plot),
between_plot_sd = sd(mean_AB, na.rm = TRUE),
between_plot_mean = mean(mean_AB, na.rm = TRUE),
.groups = "drop"
)
within_plot_summary <- pairs_complete %>%
dplyr::group_by(response, Bodentiefe) %>%
dplyr::summarise(
n_pairs = n(),
within_AB_abs_mean   = mean(abs_diff_AB, na.rm = TRUE),
within_AB_abs_median = median(abs_diff_AB, na.rm = TRUE),
within_AB_sd         = sd(abs_diff_AB, na.rm = TRUE),
.groups = "drop"
)
pseudorep_diagnostic <- within_plot_summary %>%
dplyr::left_join(
between_plot_summary,
by = c("response", "Bodentiefe")
) %>%
dplyr::mutate(
ratio_within_between = within_AB_abs_mean / between_plot_sd
)
pseudorep_diagnostic
####### BRT #########
predictors_common <- c(
"mean_temp_simple_Mar.Dec",
"Bodentiefe",
"study_site"
)
analysis_root <- analysis_root %>%
dplyr::left_join(water_level_parameters, by = c("Plot", "study_site"))
analysis_rhizome <- analysis_rhizome %>%
dplyr::left_join(water_level_parameters, by = c("Plot", "study_site"))
analysis_accumulation <- analysis_accumulation %>%
dplyr::left_join(water_level_parameters, by = c("Plot", "study_site"))
analysis_pfp <- analysis_pfp %>%
dplyr::left_join(water_level_parameters, by = c("Plot", "study_site"))
##### lmer ########
mean(basis_data$Prozentualer_Gewichtsverlust_Wurzeln, na.rm = T)
mean(basis_data$BM_accumulation_belowground_per_m2, na.rm = T)
mean(basis_data$peat_formation_potential_depth_wise, na.rm = T)
model_dec <- lmer(
mittlerer_prozentualer_Gewichtsverlust_Wurzeln_Rhizome ~ Bodentiefe +
(1 | study_site/Plot),
data = basis_data
)
summary(model_dec)
plot(basis_data$Prozentualer_Gewichtsverlust_Wurzeln~ basis_data$Bodentiefe)
qqPlot(resid(model_dec), main = "QQ-Plot: Residuen (LMM)")
plot(fitted(model_dec), resid(model_dec), xlab = "fitted", ylab = "residuals")
model_dec_roots <- lmer(
Prozentualer_Gewichtsverlust_Wurzeln ~ Bodentiefe +
(1 | study_site/Plot),
data = basis_data
)
summary(model_dec_roots)
qqPlot(resid(model_dec_roots), main = "QQ-Plot: Residuen (LMM)")
plot(fitted(model_dec_roots), resid(model_dec_roots), xlab = "fitted", ylab = "residuals")
model_dec_rhizomes <- lmer(
Prozentualer_Gewichtsverlust_Rhizome ~ Bodentiefe +
(1 | study_site/Plot),
data = basis_data
)
summary(model_dec_rhizomes)
qqPlot(resid(model_dec_rhizomes), main = "QQ-Plot: Residuen (LMM)")
plot(fitted(model_dec_rhizomes), resid(model_dec_rhizomes), xlab = "fitted", ylab = "residuals")
model_prod <- lmer(
log(BM_accumulation_belowground_per_m2) ~ Bodentiefe +
(1 | study_site/Plot),
data = basis_data
)
isSingular(model_prod, tol = 1e-5)
VarCorr(model_prod)
summary(model_prod)
qqPlot(resid(model_prod), main = "QQ-Plot: Residuen (LMM)")
plot(fitted(model_prod), resid(model_prod), xlab = "fitted", ylab = "residuals")
model_pfp <- lmer(
log(peat_formation_potential_depth_wise) ~ Bodentiefe +
(1 | study_site/Plot),
data = basis_data
)
isSingular(model_pfp, tol = 1e-5)
VarCorr(model_pfp)
summary(model_pfp)
qqPlot(resid(model_pfp), main = "QQ-Plot: Residuen (LMM)")
plot(fitted(model_pfp), resid(model_pfp), xlab = "fitted", ylab = "residuals")
####### means ######
basis_data %>%
group_by(Plot, Pseudoreplikat) %>%
summarise(
dec_combi = mean(mittlerer_prozentualer_Gewichtsverlust_Wurzeln_Rhizome, na.rm = TRUE),
.groups = "drop"
) %>%
summarise(
mean = mean(dec_combi, na.rm = TRUE),
sd   = sd(dec_combi, na.rm = TRUE),
n    = sum(!is.na(dec_combi))
)
basis_data %>%
group_by(Plot, Pseudoreplikat) %>%
summarise(
dec_combi_roots = mean(Prozentualer_Gewichtsverlust_Wurzeln, na.rm = TRUE),
.groups = "drop"
) %>%
summarise(
mean = mean(dec_combi_roots, na.rm = TRUE),
sd   = sd(dec_combi_roots, na.rm = TRUE),
n    = sum(!is.na(dec_combi_roots))
)
basis_data %>%
group_by(Plot, Pseudoreplikat) %>%
summarise(
dec_combi_rhizomes = mean(Prozentualer_Gewichtsverlust_Rhizome, na.rm = TRUE),
.groups = "drop"
) %>%
summarise(
mean = mean(dec_combi_rhizomes, na.rm = TRUE),
sd   = sd(dec_combi_rhizomes, na.rm = TRUE),
n    = sum(!is.na(dec_combi_rhizomes))
)
basis_data %>%
filter(Bodentiefe == 1) %>%
group_by(Plot, Pseudoreplikat) %>%
summarise(
dec_1_roots = mean(Prozentualer_Gewichtsverlust_Wurzeln, na.rm = TRUE),
.groups = "drop"
) %>%
summarise(
mean = mean(dec_1_roots, na.rm = TRUE),
sd   = sd(dec_1_roots, na.rm = TRUE),
n    = sum(!is.na(dec_1_roots))
)
basis_data %>%
filter(Bodentiefe == 2) %>%
group_by(Plot, Pseudoreplikat) %>%
summarise(
dec_2_roots = mean(Prozentualer_Gewichtsverlust_Wurzeln, na.rm = TRUE),
.groups = "drop"
) %>%
summarise(
mean = mean(dec_2_roots, na.rm = TRUE),
sd   = sd(dec_2_roots, na.rm = TRUE),
n    = sum(!is.na(dec_2_roots))
)
basis_data %>%
filter(!is.na(BM_accumulation_belowground_per_m2)) %>%
group_by(Plot, Pseudoreplikat) %>%
summarise(
bm_sum = sum(BM_accumulation_belowground_per_m2),
.groups = "drop"
) %>%
summarise(
mean = mean(bm_sum, na.rm = TRUE),
sd   = sd(bm_sum, na.rm = TRUE),
n    = sum(!is.na(bm_sum))
)
basis_data %>%
group_by(study_site, Plot, Pseudoreplikat, Bodentiefe) %>%
summarise(
bm_rep_depth = mean(BM_accumulation_belowground_per_m2, na.rm = TRUE),
.groups = "drop"
) %>%
group_by(study_site, Bodentiefe) %>%
summarise(
mean_bm_g_m2 = mean(bm_rep_depth, na.rm = TRUE),
sd_bm_g_m2   = sd(bm_rep_depth, na.rm = TRUE),
n            = sum(!is.na(bm_rep_depth)),
.groups      = "drop"
)
basis_data %>%
group_by(study_site, Plot, Pseudoreplikat, Bodentiefe) %>%
summarise(
bm_rep_depth = mean(BM_accumulation_belowground_per_m2, na.rm = TRUE),
.groups = "drop"
) %>%
group_by(Bodentiefe) %>%
summarise(
mean_bm_g_m2 = mean(bm_rep_depth, na.rm = TRUE),
sd_bm_g_m2   = sd(bm_rep_depth, na.rm = TRUE),
n            = sum(!is.na(bm_rep_depth)),
.groups      = "drop"
)
basis_data %>%
group_by(Plot, Pseudoreplikat) %>%
summarise(
PFP_g_m2 = mean(peat_formation_potential_depth_wise, na.rm = TRUE),
.groups = "drop"
) %>%
summarise(
mean = mean(PFP_g_m2, na.rm = TRUE),
sd   = sd(PFP_g_m2, na.rm = TRUE),
n    = sum(!is.na(PFP_g_m2))
)
basis_data %>%
group_by(Plot, Pseudoreplikat, Bodentiefe) %>%
summarise(
PFP_rep_depth = mean(peat_formation_potential_depth_wise, na.rm = TRUE),
.groups = "drop"
) %>%
group_by(Bodentiefe) %>%
summarise(
mean = mean(PFP_rep_depth, na.rm = TRUE),
sd   = sd(PFP_rep_depth, na.rm = TRUE),
n    = sum(!is.na(PFP_rep_depth)),
.groups = "drop"
)
###### ABOVEGROUND Dec und PROD BRT ###
production_above <- read.csv("production_aboveground.csv", sep = ";")
decomposition_above <- read.csv("Decomposition_aboveground.csv", sep = ";")
env
above_production_env <- inner_join(production_above, env, by = intersect(names(production_above), names(env)))
above_decomposition_env <- inner_join(decomposition_above, env, by = intersect(names(decomposition_above), names(env)))
above_production_env$study_site <- as.factor(above_production_env$study_site)
above_production_env$cover_Phragmites_.percent <- as.numeric(above_production_env$cover_Phragmites_.percent)
above_decomposition_env$Pseudoreplikat <- as.factor(above_decomposition_env$Pseudoreplikat)
above_decomposition_env$study_site <- as.factor(above_decomposition_env$study_site)
predictors_above <- c("mean_annual_water_level_m", "study_site", "mean_temp_simple_Mar.Dec" )
#BRT
above_production_clean <- above_production_env[!is.na(above_production_env[, 8]), ]
mod_prod_above <- gbm.step(data = above_production_clean, gbm.x = predictors_above, gbm.y = 8, family = "gaussian",tree.complexity = 2, learning.rate = 0.01, bag.fraction = 0.9)
##### GAM ABOVEGROUND - BELOWGROUND PRODUCTION #####
belowground_production_plot <- accumulation_response_depth %>%
dplyr::group_by(Plot, study_site) %>%
dplyr::summarise(
BM_accumulation_belowground_0_30 =
sum(BM_accumulation_belowground_per_m2, na.rm = TRUE),
.groups = "drop"
)
aboveground_production_plot <- production_above %>%
dplyr::select(
Plot, study_site, dw_reed_g
)
above_below_production <- belowground_production_plot %>%
dplyr::inner_join(
aboveground_production_plot,
by = c("Plot", "study_site")
) %>%
dplyr::filter(
!is.na(BM_accumulation_belowground_0_30),
!is.na(dw_reed_g)
)
mod_above_below_production <- mgcv::gam(
BM_accumulation_belowground_0_30 ~ s(dw_reed_g, k = 4),
data = above_below_production,
family = gaussian(),
method = "REML"
)
summary(mod_above_below_production)
AIC(mod_above_below_production)
gam.check(mod_above_below_production)
ggplot(
above_below_production,
aes(
x = dw_reed_g,
y = BM_accumulation_belowground_0_30
)
) +
geom_point(size = 2) +
geom_smooth(
method = "gam",
formula = y ~ s(x, k = 4),
method.args = list(method = "REML"),
se = TRUE
) +
labs(
x = "Aboveground biomass production (g m-2)",
y = "Belowground biomass production 0-30 cm (g m-2)"
) +
theme_bw()
##### GAM ABOVEGROUND - BELOWGROUND PRODUCTION #####
names(accumulation_response_depth)
nrow(accumulation_response_depth)
accumulation_response_depth %>%
dplyr::select(
Plot,
study_site,
Bodentiefe,
BM_accumulation_belowground_per_m2
) %>%
dplyr::arrange(Plot, Bodentiefe)
accumulation_response_depth %>%
dplyr::group_by(Plot, Bodentiefe) %>%
dplyr::summarise(
n = dplyr::n(),
mean_bm = mean(BM_accumulation_belowground_per_m2, na.rm = TRUE),
.groups = "drop"
)
##### GAM ABOVEGROUND - BELOWGROUND PRODUCTION #####
unique(basis_data$Bodentiefe)
basis_data %>%
dplyr::group_by(Bodentiefe) %>%
dplyr::summarise(
mean_bm = mean(BM_accumulation_belowground_per_m2, na.rm = TRUE),
n = sum(!is.na(BM_accumulation_belowground_per_m2)),
.groups = "drop"
)
##### GAM ABOVEGROUND - BELOWGROUND PRODUCTION #####
read.csv("accumulation_belowground_per_m2_0-30.csv")
##### GAM ABOVEGROUND - BELOWGROUND PRODUCTION #####
prod_below <- read.csv("accumulation_belowground_per_m2_0-30.csv")
prod_below
belowground_production_depth <- prod_below %>%
dplyr::filter(!is.na(BM_accumulation_belowground_per_m2)) %>%
dplyr::group_by(Plot, Bodentiefe) %>%
dplyr::summarise(
BM_accumulation_belowground_per_m2 =
mean(BM_accumulation_belowground_per_m2, na.rm = TRUE),
.groups = "drop"
)
belowground_production_plot <- belowground_production_depth %>%
dplyr::group_by(Plot) %>%
dplyr::summarise(
BM_accumulation_belowground_0_30 =
sum(BM_accumulation_belowground_per_m2, na.rm = TRUE),
n_depths = dplyr::n(),
.groups = "drop"
) %>%
dplyr::filter(n_depths == 3)
aboveground_production_plot <- production_above %>%
dplyr::group_by(Plot, study_site) %>%
dplyr::summarise(
dw_reed_g = mean(dw_reed_g, na.rm = TRUE),
.groups = "drop"
) %>%
dplyr::filter(!is.na(dw_reed_g))
above_below_production <- belowground_production_plot %>%
dplyr::inner_join(
aboveground_production_plot,
by = "Plot"
)
mod_above_below_production <- mgcv::gam(
BM_accumulation_belowground_0_30 ~ s(dw_reed_g, k = 4),
data = above_below_production,
family = gaussian(),
method = "REML"
)
summary(mod_above_below_production)
AIC(mod_above_below_production)
gam.check(mod_above_below_production)
ggplot(
above_below_production,
aes(
x = dw_reed_g,
y = BM_accumulation_belowground_0_30
)
) +
geom_point(size = 2) +
geom_smooth(
method = "gam",
formula = y ~ s(x, k = 4),
method.args = list(method = "REML"),
se = TRUE
) +
labs(
x = "Aboveground biomass production (g m-2)",
y = "Belowground biomass production 0-30 cm (g m-2)"
) +
theme_bw()
